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1 1. Models of hybrid incompatibility 

2 

3 In order to investigate changes in allele frequencies at hybrid incompatibility loci, we 

4 need to define explicit models of selection on each two-locus genotype. Different genotypes at 

5 hybrid incompatibility loci will experience different strengths of selection. This can be a 

6 consequence of dominance at each locus (h), the number of epistatic interactions, the order in 

7 which mutations occurred, and the relative fitness of the parental and ancestral genotypes. 

8 Our simulations are based on fitness matrices representing adaptive BDM 

9 incompatibilities (Fig. SI) and coevolving hybrid incompatibilities (Fig. S2). Elsewhere, we 

10 discuss neutral BDM incompatibilities and coevolving hybrid incompatibilities that will not 

1 1 generate isolation in the absence of strong genetic drift (Fig. S7). For simplicity, in most cases 

12 we simulate coevolving incompatibilities (Fig. S2) and assume that selection on different 

13 epistatic interactions is symmetrical (si=S2) and that all alleles are codominant (h=0.5). However, 

14 we relax these assumptions in Supporting Information 3C and D. Though we focus on 

15 coevolving incompatibilities (Fig. S2) our results also generally apply to adaptive BDM 

16 incompatibilities (Fig. SI) which have a similar fitness matrix. 

17 In Figures SI A and S2A, we show a single mutational path to each hybrid 

18 incompatibility type; however, in both cases, several possible mutational paths exist. Because the 

19 lineage and order in which mutations occur is expected to be random (1), we randomly assign a 

20 fitness matrix to each incompatibility in simulations (note: these matrices are identical at h=0.5 

21 and si^). 

22 Using these fitness matrices and Karlin's two-locus selection model (see Methods, (2)), 

23 we can determine the expected trajectory of incompatibility pairs in hybrid populations under 



1 different parameters (Fig. S3). However, some of these dynamics differ in finite populations (see 

2 next section). 
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A) Two-locus BDM model 
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5 

6 Fig. SI. Evolution of two-locus BDM incompatibility and fitness of hybrid genotypes. 

7 (A) One of two possible mutational paths to the development of a two-locus BDM incompatibility 

8 (Not shown is the case where both mutations occur on one lineage). These incompatibilities can arise 

9 as the result of neutral fixation (Wp ar entai=w ancestr ai) or as the result of adaptive evolution 

10 (Wparenta^Wancestrai)- (B) Potential selection patterns on hybrid genotypes between the two parentals 

1 1 (assuming Wparenta^Wancestrai); genotypes corresponding to selection coefficients s\ and are indicated 

12 in blue and red respectively. For BDM incompatibilities, si and 52 can be asymmetric, and in neutral 

13 BDM incompatibilities either si or s 2 will equal zero. (C) Fitness of hybrid individuals with each 

14 genotype will depend the intensity of selection (s\, si) and dominance Qia, h B ) at the two loci. We 

15 assume, for simplicity that the fitness advantage of all derived genotypes (here, xB and Ax) is equal. 
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A) Two-locus co-evolutionary model 
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B) Possible epistatic interactions 
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Fig. S2. Evolution of two-locus coevolving incompatibility and fitness of hybrid genotypes. (A) 

One of two possible mutational paths to the development of a two-locus coevolved incompatibility. 
Not shown is the case where B 2 precedes A (see Supplement 4). (B) Potential epistatic interactions 
among hybrid genotypes. Incompatibilities corresponding to S\ and S2 are indicated in blue and red, 
respectively. (C) Fitness of hybrid individuals with each genotype will depend the intensity of 
selection (si, si) and dominance (h A , h B ) at the two loci. 
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A) Neutral BDM model 



B) Coevolution model: codominant symmetrical 
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C) Coevolution model: codominant asymmetrical 
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E) Coevolution model: h A = h B1 = 1 
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D) Coevolution model: h A # h E 
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F) Coevolution model: h A = h B1 = 0 



00 

d 



CD 

o 



d 



o 
d 




1 Generations Generations 

2 Fig. S3. Fixation of two-locus hybrid incompatibilities under deterministic selection (the Karlin 

3 model). The expected parent 1 -derived allele trajectories for two unlinked hybrid incompatibilities 

4 under the Karlin model depend on starting admixture proportions (/=0.3-0.7 shown here), dominance 

5 parameters (h), and the intensity of selection (i.e. Si, si, see Figs. SI and S2). The solid line tracks 

6 ancestry at locus 1 and the dashed line shows ancestry at locus 2. (A) Neutral BDM incompatibility 

7 pairs do not fix iff 0.5; at /= 0.5 they fix for a hybrid genotype pair that is not incompatible with 

8 either parental species (see Supporting Information 4). When incompatibilities are codominant (B, C), 

9 the two incompatibility loci fix deterministically for the major parent. At certain values of h (D, E, F), 
10 fixation is less dependent on initial admixture proportions. 

11 
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1 2. Karlin model with genetic drift (multinomial sampling) and comparison to 

2 population simulations. 

3 The Karlin model (ref. (2); see Methods) allows us to predict the patterns of fixation for 

4 different incompatibility types with different parameters (Fig. S3) but is not realistic because 

5 even large populations will have some genetic drift, which is not accounted for in the model. To 

6 implement drift, we added multinomial sampling of N individuals at each generation (Fig. S4). 

7 The major difference caused by implementing drift is that fixation patterns do not always follow 

8 the predictions of the Karlin model. This can be seen most clearly in Fig. S4C, where some 

9 fixation trajectories differ from expectations under the Karlin model. 

10 The Karlin equations do not allow us to model the fates of more than one incompatibility 

1 1 pair. We thus wrote an explicit population simulation code (admix ' em) described in the 

12 Methods. To validate this code, we compared results from simulations with admix ' em to 

13 simulations of Karlin' s model with drift implemented by multinomial sampling. We simulated a 

14 single coevolving hybrid incompatibility (h=0.5, 51=52=0.1) with population size N=1,000 and 

15 10,000. In both cases the results from admix ' em and Karlin's model with drift are nearly 

16 identical (Table SI). 
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A) Karlin model B) Karlin model with N=1 0,000 C) Karlin model with N=1 ,000 




0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 

j Generations Generations Generations 

2 Fig. S4. Fixation of hybrid incompatibilities under the Karlin model with genetic drift. The 

3 expected patterns of fixation for a coevolving hybrid incompatibility (Fig. S2) under Karlin's model 

4 depend on starting admixture proportions (/=0.3-0.7 shown here), dominance parameters (h), and 

5 asymmetry in selection (s\ ^ S2). The solid line shows ancestry at locus 1 of an incompatibility and 

6 the dashed line shows ancestry at locus 2 of an incompatibility. (A) Parent 1 allele trajectories 

7 predicted by the Karlin model for a given set of parameters. (B) Results for the same parameters 

8 incorporating multinomial sampling of 10,000 individuals at each generation. (C) Results for the 

9 same parameters incorporating multinomial sampling of 1,000 individuals at each generation. Patterns 

10 of fixation depend less on initial admixture proportions as drift increases. The equilibrium at f=0.5 in 

1 1 the selection-only model is unstable in the presence of drift. 
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1 Table SI. Comparison of rates of fixation based on the Karlin model with genetic drift and the 

2 admix ' em simulation program at different population sizes. 





Karlin model with multinomial 


admix ' 


em 




sampling 








Diploid 


Percent fixing parent 


Average time 


Percent fixing 


Average time 


population size 


1 genotype ± SE 


to fixation 


parent 1 genotype 


to fixation 


(N) 




± SD 


±SE 


± SD 


1,000 


49.8 ±2 


152 ±42 


50.2 ±2 


154 ±40 


10,000 


50.6 ±2 


201 ± 44 


51.0 ±2 


193 ± 32 


Note - One hybrid incompatibility pair (Fig. S2), si=S2=0.l, 


f=0.5, h=0.5 for 500 replicate 



4 simulations. 
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1 3. Evaluation of the hybrid population reproductive isolation model under a 

2 range of parameters 

3 

4 Simulations presented in the main text demonstrate that with two hybrid incompatibility 

5 pairs, reproductive isolation frequently evolves (47±2%) due to fixation of hybrid 

6 incompatibility loci for both parents. This mechanism of hybrid isolation is most likely to be 

7 biologically relevant if it occurs under a broad range of parameters. We performed simulations of 

8 two coevolving hybrid incompatibility pairs (Fig. S2) varying selection (s), dominance (h), 

9 diploid population size (N), admixture proportions if), and degree of asymmetry in selection 

10 {s\£s2). All results are based on 500 replicate simulations. We consider hybrid populations 

1 1 "isolated" from both parents if they fix for at least one incompatibility pair from each parent; the 

12 strength of isolation will always depend on the strength of selection on each incompatibility and 

13 the number of incompatibilities. Together, these simulations demonstrate that hybrid 

14 reproductive isolation due to fixation of genetic incompatibilities can occur under a range of 

15 parameters. 
16 

17 A. The strength of selection on hybrids. 

18 Variation in the strength of selection on hybrids can influence the likelihood of hybrid 

19 isolation under our model. We performed simulations with two coevolving hybrid 

20 incompatibility pairs with different strengths of selection (Table S2, Fig. S5). We find that the 

21 probability of isolation decreases as total selection on hybrids increases. Lower probability of 

22 isolation at higher levels of selection on hybrids is the result of hybrid populations becoming 

23 dominated by parental individuals. 
24 
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Fig. S5. Probability of isolation with increasing selection on hybrids. 

With increasing selection on Fi hybrids between the parental species, the 
probability that hybrid populations will develop reproductive isolation 
from both parents decreases. This result is less pronounced with a greater 
number of hybrid incompatibility (HI) pairs. Error bars show two standard 
deviations. Simulation parameters were h=0.5, Si=Si,f=0.5, and N=1,000. 



Table S2. The effect of increasing selection on hybrids on the probability of and time to 
isolation. 



Fitness of Fl 
hybrid 



Percent isolating ± SE 



Average time to 
isolation ± SD 



15 
16- 



0.9 
0.8 
0.7 
0.5 



47 ±2 
36 ±2 
27 ±2 
14 ±2 



150 ±48 
78 ±21 
55 ± 14 
53 ±9 



17 
18 
19 
20 
IX. 



Note - Two hybrid incompatibility pairs (Fig. S2), si=s 2 , N=1000,/=0.5, h=0.5 for 500 
replicate simulations. 
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1 B. Population size 

2 A number of models of speciation rely on small population sizes to achieve isolation from 

3 parental species (e.g. see (3, 4)). One feature of our model is that it does not rely on small 

4 population sizes; in large populations incompatibility pairs fix by selection, whereas in small 

5 populations they fix by chance. As a result, the proportion of hybrid populations evolving 

6 isolation from both parents is not strongly dependent on population size (Table S3). This 

7 suggests that our model is applicable to a range of natural hybrid populations. 
8 

9 Table S3. The effect of population size on the probability of 

10 and time to isolation. 



Diploid 


Percent 


Average time to 


population size 


isolating 


isolation 


(N) 


± SE 


± SD 


100 


40 ±2 


75 ±34 


1000 


47 ±2 


150 ±48 


10000 


44 ±2 


190 ± 43 



1 1 Note - Two hybrid incompatibility pairs (Fig. S2), si=s 2 =0. 1 , 

12 ,/=0.5, /z=0.5 for 500 replicate simulations. 
13 

14 
15 

16 C. Initial admixture proportions 
17 

18 To determine how changing admixture proportions influences the number of hybrid 

19 populations evolving isolation, we varied starting admixture proportions fromy=0.5-0.7. As 

20 expected, the probability of isolation from both parents decreased with increasing ancestry skew 

21 (Fig. S6). Initial admixture proportions close/=0.5 are most likely to result in isolation from both 

22 parental species. However, even in skewed populations, isolation from both parents was possible 

23 when assuming variation in dominance among alleles (for e.g. h =Uniform(0,l), Fig. S6) or at 

24 lower population sizes (N=100, Fig. S6). 
25 
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Fig. S6. The effect of initial admixture proportion on the 
probability of isolation. Proportion of hybrid populations 
developing isolation from both parents as a function admixture 
proportions, dominance (h) and population size (two 
incompatibility pairs, sl=s2=0.1). Isolation occurs most 
frequently at equal admixture proportions, but can occur in 
ancestry-skewed populations, especially if the populations are 
small or there is variation in dominance. Error bars show two 
standard deviations. 



D. Asymmetric selection on hybrid incompatibilities 

For simplicity, in most cases we simulate symmetric selection on hybrid 
incompatibilities. To introduce asymmetric selection on different hybrid genotypes we use 
different selection coefficients for distinct interactions in the coevolution model, S\ and S2 (Fig. 
S2). When selection is asymmetrical (i.e. S\ * s?), the Karlin model has similar properties to S\ = 
S2 as long as Si and S2 are both greater than 0 (Fig. S3). To investigate this in finite populations, 
we randomly assigned 5=0.1 to s\ or S2 independently for two coevolving incompatibility pairs 
and set the value for the other selection coefficient to 0.2*5, 0.3*5, 0.4*5, or 0.5*5. In these 
simulations, reciprocal isolation still developed frequently (Table S4). 



1? 



1 The adaptive BDM model (Fig. SI) might result in even more extreme asymmetry in 

2 selection on hybrid genotypes (5-7). To approximate this scenario, we simulated adaptive 

3 BDMIs (Fig. SI) and randomly assigned 5=0.5 to S\ or s 2 independently for two hybrid 

4 incompatibility pairs, and assigned 5=0.02 to the other selection coefficient. In these simulations 

5 39 ± 2% of populations became isolated from both parental species. In populations fixing 

6 incompatibility pairs from both parents, the strength of isolation from each parent species was 

7 similar. On average, the absolute fitness difference between backcross types (hybrid x parent 1 

8 and hybrid x parent 2) was only 0.003. The average fitness of a hybrid x parent cross after hybrid 

9 populations fixed for parental incompatibilities was 0.87 (compared to 0.76 for an Fl cross 
10 between the two parent species). 

11 

12 Table S4. Effects of asymmetry in selection 

13 on the probability of isolation. 



Asymmetry 


Percent 


Average time 


of selection 


isolating 


to isolation ± 




± SE 


SD 


sa=sb 


47 ±2 


150 ±48 


s a =0.5*sb 


44 ±2 


211 ±64 


sa=0.4*s b 


46 ±2 


232 ± 73 


s A =0.3*s B 


44 ±2 


290 ± 86 


s a =0.2*s b 


36 ±2 


319 ±85 


s a =0*sb 


0.2 ± 0.2 


497 (NA) 



14 Note - Two hybrid incompatibility pairs 

15 (Fig. SI, S2), s A =0.1, N=1000,y=0.5, h=0.5 

16 for 500 replicate simulations. 
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1 E. Dominance 
2 

3 In the above simulations, we model codominant hybrid incompatibilities but varying 

4 dominance (h) does not have major effects on our conclusions. To demonstrate this, we 

5 simulated two scenarios: 1) we drew h from a uniform distribution, requiring h at each locus to 

6 add to 1 and 2) we randomly assigned h discrete values of 0 (recessive), 0.5 (co-dominant), or 1 

7 (dominant), again requiring h to add to 1 at each locus. Both of these scenarios result in a similar 

8 probability of isolation from parental species (Table S5). 
9 

10 Table S5. The effect of variation in dominance. 



Dominance 


Percent 


Average time to 




isolating 


isolation ± SD 




± SE 




A = 0.5 


47 ±2 


150 ±48 


/z = Uniform(0,l) 


48 ±2 


120 ±36 


h G (0, 0.5, 1) 


44 ±2 


118 ± 37 



1 1 Note - Two hybrid incompatibility pairs (Fig. S2), 

12 5i=52=0.1, N=1000,/=0.5, for 500 replicate simulations. 
13 
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1 4. Incompatibilities that do not frequently result in reproductive isolation in 

2 the absence of strong drift 

3 

4 Several types of genetic incompatibilities, such as neutral BDM incompatibilities (or 

5 other scenarios that generate identical fitness matrices, Fig. S7), are not predicted to contribute to 

6 reproductive isolation between hybrids and parentals by incompatibility selection (Fig. S7). This 

7 is because hybrid incompatibility pairs either do not fix, or fix for a genotype that is compatible 

8 with both parents (based on the Karlin model, Fig. S3A, Fig. S7). Simulations of 

9 incompatibilities with these fitness matrices confirm this prediction; the most frequently fixed 

10 genotype is not incompatible with either parent (54±2% of simulations) and in many simulations 

1 1 (45±2%) hybrid populations did not fix for a particular parental genotype at incompatibility loci 

12 within 2,000 generations (simulation parameters: N=1000, 5=0.1,^=0.5, h=0.5). 

13 Nonetheless, with strong genetic drift, genotypes incompatible with one of the parents 

14 can fix by chance, contributing to isolation, even under a neutral BDMI model. This can be seen 

15 in Fig. S8, which shows the proportion of hybrid populations isolated from both parents as a 

16 function of population size. As population size decreases and drift increases, the proportion of 

17 hybrid populations isolating from both parents increases. However, the time to isolation is longer 

18 than in cases where fixation is driven by deterministic selection (by drift N=100: 215±105 

19 generations, N=500: 1267±506 generations, N=1000: 1312±278 generations, compared to Table 

20 S3). 
21 
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A) Neutral BDM incompatibilities do not result in isolation 
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B) Co-evolving incompatibilities that do not result in isolation 
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2 Fig. S7. Hybrid incompatibility models that do not frequently result in reproductive isolation 

3 between hybrid and parental populations in the absence of drift. (A) When hybrid populations form 

4 at equal admixture proportions, the deterministic model predicts that neutral BDM incompatibilities will 

5 fix for the ancestral genotype in a two-lineage model (left) and a genotype that is compatible with both 

6 species in a one-lineage model (right). (B) In a coevolution scenario, certain mutation orders result in an 

7 identical fitness matrix to A and thus do not result in reproductive isolation in hybrid populations. In all 

8 cases depicted, mutations in lineage 1 could occur in lineage 2 and vice versa but the expected effects on 

9 isolation from parental species do not change. 
10 
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3 Fig. S8. The effect of population size on the 

4 probability of isolation due to neutral BDM 

5 incompatibilities. As drift increases, the proportion 

6 of hybrid populations isolated from parentals by 

7 fixation of neutral BDM incompatibilities increases. 

8 However, this process does not occur as rapidly as 

9 deterministic selection on other types of hybrid 

10 incompatibilities. Simulation parameters: two 

1 1 neutral BDMI pairs (Fig. S7), 5=0. 1 ,f=0.5, h=0.5 

12 for 500 replicate simulations. 
13 
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1 5. Validation of the model under a range of genetic architectures 

2 

3 Because species are likely to be differentiated by a number of incompatibilities even in the 

4 early stages of speciation (8), we investigated how the number of incompatibility pairs and their 

5 genetic architecture affects the probability of evolving reproductive isolation. 
6 

7 A. Number of incompatibility pairs 

8 Theoretically, increasing the number of incompatibility pairs should increase the 

9 probability that hybrid populations will become isolated from both parents (by at least one 

10 incompatibility). However, increasing the number of incompatibility pairs will reduce hybrid 

1 1 fitness and thus may also increase the probability that hybrids will be outcompeted by parental 

12 species. Simulating 3-6 unlinked coevolving incompatibility pairs (/=0.5, h=Q.5, 51=52=0.1, 

13 N=1000) increased the probability of isolation as a function of the number of pairs (Fig. 4, 58- 

14 96%). This result is notable because the probability of reciprocal isolation increases even with an 

15 increase in total selection on hybrids (see also Fig. S5). 
16 

17 B. Probability of reciprocal isolation in the presence of neutral BDM incompatibilities 

1 8 Though there is little data on fitness matrices for specific hybrid incompatibilities 

19 (reviewed in 9), neutral BDM incompatibilities may be common (1). Because of this, we 

20 determined the frequency of isolation between hybrids and parents in simulations including 

21 neutral BDM incompatibilities, where either s\ or s 2 equals zero. We simulate two coevolving 

22 hybrid incompatibility pairs as described above but also include one linked or unlinked neutral 

23 BDM incompatibility pair. We find that the presence of BDM incompatibilities does not 

24 significantly decrease the probability of isolation between hybrid and parental populations 
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1 (unlinked: 45 ± 2%, 1 cM between coevolving incompatibility and neutral BDM as in Figure 

2 S9B: 48 ± 2%, compared to 47 ± 2% with no neutral BDMIs). Interestingly, linkage between a 

3 neutral BDM incompatibility and a coevolving incompatibility can result in more frequent 

4 fixation of one of the parental genotypes at BDMIs (21± 2%, compared to 0.2% without 

5 linkage). In such cases, neutral BDMIs will contribute to isolation between hybrids and parental 

6 species. 
7 

8 C. Linkage between incompatibility pairs 

9 Linkage between sites involved in hybrid incompatibilities could influence the frequency 

10 of isolation because ancestry at linked sites will be correlated. In order to investigate this we 

1 1 simulated two coevolving hybrid incompatibility pairs and varied genetic distance between loci 

12 (50 cM, 10 cM, and 1 cM, Fig. S9). Some scenarios of physical linkage did not have a strong 

13 effect on the probability of isolation, while others reduced the probability of isolation (Table S6). 
14 
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A) Scenario 1: Linkage between members of the same incompatibility pair 



Species 1 



X 



Species 2 



TT 



TV 



]- pair 1 



pair 2 



B) Scenario 2: Linkage between members of different incompatibility pairs 
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Fig. S9. Linkage between incompatibility pairs. Linkage between 
incompatibility pairs can change the probability of hybrid populations 
evolving reproductive isolation (Table S6). (A) In scenario 1, linkage 
between loci in the same incompatibility pair does not influence the 
frequency of hybrid populations evolving reproductive isolation. 
(B) In linkage scenario 2, linkage between loci in different incompatibility 
pairs significantly decreases the frequency at which hybrid populations 
evolve reproductive isolation. The probability of recombination between 
two sites is indicated as r. 



Table S6. Effect of linkage between hybrid incompatibilities on 
the probability of and time to isolation. 


Linkage 


Genetic 


Percent 


Average time 


scenario 


distance 


isolating ± SE 


to isolation ± 


(Fig. S9) 






SD 


Scenario 1 


50 cM 


47 ±2 


133 ±35 


Scenario 2 


50 cM 


43 ±2 


137 ±51 


Scenario 1 


10 cM 


28 ±2 


130 ±31 


Scenario 2 


10 cM 


43 ±2 


137 ±35 


Scenario 1 


1 cM 


6± 1 


157 ±23 


Scenario 2 


1 cM 


47 ±2 


122 ±28 



Note - A visual representation of the linkage scenarios can be found 
in Fig. S9. Simulation parameters: Two hybrid incompatibility pairs 
(Fig. S2), 5i=52=0.1. N=1000,y=0.5, h=0.5 for 500 replicate " 
simulations. 
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1 D. More complex incompatibility architectures 
2 

3 We have focused on very simple models that show that hybrid reproductive isolation can 

4 evolve in principle. However, empirical studies suggest that even closely related species are 

5 sometimes differentiated by dozens of incompatibilities that vary in their fitness effects (10-14). 

6 To investigate whether hybrid reproductive isolation evolves with more complex incompatibility 

7 architectures, we simulated a large number of incompatibilities with variation in selection, 

8 dominance, and asymmetry. We simulate 20 pairs of hybrid incompatibilities with a randomly 

9 chosen genomic position (chromosome 1-20, position 1-25000000 bp). Dominance is drawn 

10 from a uniform distribution (0-1) and selection coefficients (si and S2) are drawn from an 

1 1 exponential distribution with a mean of 5=0.05. Because S\ and S2 are determined independently 

12 this introduces asymmetry in selection, and if si or s 2 is close to 0 (i.e. on the order of 1/N), this 

13 generates an effectively neutral BDM incompatibility. 

14 In these simulations, hybrid populations had a high probability of developing isolation 

15 from both parental species (95±1%). On average, backcrosses between parentals and individuals 

16 from the hybrid population had a fitness of 0.79±0.1 1 after 500 generations of selection on 

17 hybrid populations, compared to an average of 0.61±0.07 for Fi hybrids between the two 

18 parental species (Figs. 3, S10). The average fitness of hybrids mating with other hybrids in the 

1 9 population was 0 . 97±0 . 02 . 
20 
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A) B) C) 
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2 Fig. S10. Mating with parents reduces fitness of allopatrically evolved hybrid populations. (A) 

3 Change in average hybrid population fitness over time in a simulation of 20 incompatibility pairs with 

4 dominance and selection coefficients drawn from an exponential distribution (see Supporting information 

5 5D). (B) The same hybrid population with a one generation burst of migrants from parent 1 (4Nmi=400) 

6 at generation 300. (C) The same hybrid population with a one generation burst of migrants from parent 2 

7 (4Nm 2 =400) at generation 300. Notably, hybrid populations have lower average fitness after gene flow 

8 with either parent, but recover rapidly. 
9 

10 
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1 6. Validation of the model under a range of demographic scenarios 

2 

3 A. Isolation with migration 

4 Hybrid populations in nature frequently experience ongoing gene flow with parental 

5 species. Given this, it is important to determine the impact of gene flow on how frequently 

6 hybrid populations evolve reproductive isolation. We simulated two migration scenarios, one 

7 with continuous migration from both parental populations (4Nm = 8 - 40) and one with a single 

8 pulse of migration from either parent (migrating parent determined randomly) at 25, 50, and 75 

9 generations in independent simulations respectively. Because populations with ongoing 

10 migration from both parents never completely fix for a particular genotype, we treated cases in 

1 1 which 95% of the hybrid population had a particular genotype as fixed. As expected, higher 

12 migration rates lowered the proportion of hybrid populations evolving isolation from both 

13 parents (Table S7, compared to 47± 2% with no migration). However, our results demonstrate 

14 that isolation still occurs frequently with fairly high levels of gene flow from the parental 

15 species. Notably, even before fixation, pulses of migration from parental populations lower 

16 hybrid population fitness (Fig. 3, Fig. S9). 
17 

18 B. Migration in a stepping stone hybrid zone 

19 We use the simplest hybrid zone structure in our simulations (Fig. 1A of the manuscript), 

20 but hybrid zones often exist as a series of populations between the parental species. To simulate 

21 this, we modeled three hybrid populations between the parental populations with bidirectional 

22 migration between adjacent populations (Fig. SI 1, Table S7). In this scenario, the hybrid 

23 population at the center of the hybrid zone (H2, Fig. SI 1) experiences less direct gene flow from 
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1 the parental populations, and should evolve reproductive isolation more frequently at the same 

2 migration rate. As expected, we found that the hybrid population at the center of the hybrid zone 

3 developed isolation from parental species more frequently (Table S7), but that hybrid 

4 populations closer to the parentals were less likely to develop isolation (4Nm=S: 26±2 - 30±2%, 

5 4Nm=20: 1.6±0.6 - 2.0±0.6%, 4Nm=40: 0). This suggests that more realistic hybrid zone 

6 structures could actually increase the likelihood that some hybrid populations will evolve 

7 reproductive isolation by buffering effective gene flow from parental species. 
8 
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Fig. Sll. Hybrid zone structure used in simulations. Model of hybrid zone structure used in 
simulations of complex hybrid zone structures (see Supporting Information 6B). This structure of a 
gradient of hybrid populations with ongoing gene flow from parental and other hybrid populations is 
similar to many naturally occurring hybrid populations. 
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1 

2 Table S7. The effect of on-going and bursts of parental immigration 

3 on the probability of and time to isolation. 

Migration scenario Hybrid zone Percent Average time to 

structure isolating isolation ± SD 



± SE 



Continuous (4Nm=S) 


Fig. 


1A 


J 1 x z 


1 id 




4Q 


Continuous (4Nm=20) 


Fig. 


1A 


27 ±2 


247 


± 


81 


Continuous (4/Vm=40) 


Fig. 


1A 


0.6 ±0.3 


384 


± 


86 


Burst - generation 25 (4iVm=400) 


Fig. 


1A 


36 ±2 


153 


± 


48 


Burst - generation 50 (47Vra=400) 


Fig. 


1A 


38 ±2 


151 


± 


44 


Burst - generation 75 (4iVm=400) 


Fig. 


1A 


38 ±2 


152 


± 


47 


Continuous (4Nm=S) 


Fig. 


Sll 


47 ±2 


153 


± 


44 


Continuous (4^=20) 


Fig. 


Sll 


32 ±2 


233 


± 


79 


Continuous (4Nm=40) 


Fig. 


Sll 


1 ±0.4 


371 


± 


130 


Note - Two hybrid incompatibility pairs (Fig. S2), sr 


=52=0.1, N= 1000, 









5 7=0.5, h=0.5 for 500 replicate simulations. For the "Burst" models, all 

6 migration occurs in one generation. 
7 

8 C. Reciprocal isolation of multiple hybrid populations 

9 One interesting prediction of our model is that independent hybrid populations that form 

10 from the same parental species could evolve reproductive isolation from each other, if they fix 

1 1 for different parental alleles at incompatibility pairs. To investigate this, we simulated two hybrid 

12 populations formed from the same parental species (Fig. S12). In cases where both hybrid 

13 populations in a simulation evolved isolation from the parents, we asked how frequently these 

14 populations were isolated from each other. When there was no gene flow between the hybrid 

15 populations, isolation evolved frequently between these populations (2 HI: 52 ± 2% 3HI: 80 ± 

16 2%, Table S8). However, when hybrid populations were connected by gene flow, this outcome 

17 was significantly less likely (Table S8). These simulations suggest that hybridization between the 

18 same parental species could generate multiple, reproductively isolated hybrid populations, but 

19 that this outcome is more likely when gene flow is lower. 
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1 




2 

3 Fig. S12. Structure of multiple hybrid populations. Hybrid zone structure used in simulations of 

4 reciprocal hybrid isolation (Supporting Information 6B). 
5 

6 Table S8. Independently formed hybrid populations can evolve reproductive 

7 isolation from each other. 



Number of 


Migration 


Percent reciprocally 


incompatibility pairs 


rate 


isolated ± SE 


2 


4Nm~- 


=0 


52 ± 5 


2 


4Nm~- 


=8 


18 ±4 


2 


4Nm= 


=20 


0 


3 


4Nm~- 


=0 


80 ±4 


3 


4Nm~- 


=8 


66 ± 5 


3 


4Nm= 


=20 


7±3 



8 Note - Two or three hybrid incompatibility pairs (Fig. S2), 

9 5i=52=0.1, N=1000,/=0.5. Simulations conducted until 100 

10 replicates were generated in which both hybrid populations 

1 1 were isolated from parentals. 
12 

13 

14 

15 
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1 D. Non-random mating of parents in hybrid zones 

2 In our simulations, mating is random, and pure parental individuals become less common 

3 as the number of generations since initial hybridization increases. However, if parentals mate 

4 assortatively, larger numbers will exist in the hybrid population, increasing the risk of out- 

5 competition of hybrids by parents. To investigate this, we implemented a mating rule in which 

6 migrating parentals mate exclusively with conspecifics. If a conspecific is not sampled in 25, 50, 

7 or 75 attempts in separate simulations respectively, the parental individual does not mate. With 

8 increasing parental effort to mate with a conspecific, hybrid populations evolved isolation less 

9 frequently (Table S9). This suggests that hybrid populations are more likely to be outcompeted if 
10 parentals have strong preferences against mating with hybrids. 

11 

12 Table S9. Parental preferences for conspecifics reduce the frequency of 

13 hybrid reproductive isolation. 



Parental 


Percent 


Average time to 


preference 


isolatin 


g±SE 


isolation ± SD 


None 


47± 


2% 


150 ±48 


25 attempts 


42 ± 


2% 


322 ± 85 


50 attempts 


33 ± 


2% 


384 ± 79 


75 attempts 


34 ± 


2% 


328 ± 74 



14 Note - Two hybrid incompatibility pairs (Fig. S2), 4Nmi=4Nm2=S, 

15 5i=52=0.1, N=1000,/=0.5, h=0.5 for 500 replicate simulations. 
16 

17 
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